A significant quantitative trait locus on chromosome Z and its impact on egg production traits in seven maternal lines of meat-type chicken

Background Egg production is economically important in the meat-type chicken industry. To better understand the molecular genetic mechanism of egg production in meat-type chicken, genetic parameter estimation, genome-wide association analyses combined with meta-analyses, Bayesian analyses, and selective sweep analyses were performed to screen single nucleotide polymorphisms (SNPs) and other genetic loci that were significantly associated with egg number traits in 11,279 chickens from seven material lines. Results Yellow-feathered meat-type chickens laid 115 eggs at 43 weeks of age and white-feathered chickens laid 143 eggs at 60 weeks of age, with heritability ranging from 0.034–0.258. Based on meta-analyses and selective sweep analyses, one region (10.81–13.05 Mb) on chromosome Z was associated with egg number in all lines. Further analyses using the W2 line was also associated with the same region, and 29 SNPs were identified that significantly affected estimation of breeding value of egg numbers. The 29 SNPs were identified as having a significant effect on the egg number EBV in 3194 birds in line W2. There are 36 genes in the region, with glial cell derived neurotrophic factor, DAB adaptor protein 2, protein kinase AMP-activated catalytic subunit alpha 1, NAD kinase 2, mitochondrial, WD repeat domain 70, leukemia inhibitory factor receptor alpha, complement C6, and complement C7 identified as being potentially affecting to egg number. In addition, three SNPs (rs318154184, rs13769886, and rs313325646) associated with egg number were located on or near the prolactin receptor gene. Conclusion Our study used genomic information from different chicken lines and populations to identify a genomic region (spanning 2.24 Mb) associated with egg number. Nine genes and 29 SNPs were identified as the most likely candidate genes and variations for egg production. These results contribute to the identification of candidate genes and variants for egg traits in poultry. Supplementary Information The online version contains supplementary material available at 10.1186/s40104-022-00744-w.


Introduction
Throughout the twentieth century, specialized commercial chicken populations were established for egg laying and meat production [1]. Currently, layers display high levels of egg production and have been developed by Open Access † Jiqiang Ding and Fan Ying these authors have contributed equally to this work and share first authorship.  13:96 intensive artificial selection in the egg-laying industry to yield more than 320 eggs during a 52-week laying period [2,3]. Meat-type chickens have received substantial attention for their meat production, and to avoid selection for both growth traits and egg production traits in the same bird [1]. For meat-type chickens, the characteristics of body weight, growth rate, and feed conversion rate have been widely studied [4]. Although significant genetic changes may have occurred during selective breeding practices, the relationships between genotypes and phenotypes are not well characterized [2]. Many studies have shown that the heritability of egg number production. For example, such heritability was 0.13-0.45 and 0.10-0.40 for egg number production in two commercial lines of White Leghorn hens [5]. In the past three decades, many studies have explored egg production at the genetic level. Recent studies of reproductive hormones have identified 31 candidate genes significantly associated with egg production, and there are 64 new candidate genes and 108 SNPs associated with oviposition performance based on genome sequencing analyses [6]. Over 440 quantitative trait loci (QTL) on the Gallus gallus chromosome (GGA) have been reported to be associated with egg number in chickens [7]. With the deeply sequencing and the mining of increasing SNPs, genome-wide association studies (GWAS) and selective sweep analyses have become powerful tools for detecting QTLs in many breeds of livestock and poultry [8].
The growth and feed conversion traits of meat-type chickens have been widely reported, however, egg-laying traits have been less studied. For layers, multiple SNPs and QTLs have been found [6], but due to the genetic differences between layers and meat-type chickens, the application of those findings in meat-type chickens has been limited. Mapping of QTLs and SNPs with potential impact on egg production in meat-type chickens to provide more reference information for subsequent genome selection. Our experimental design involved GWAS and selective sweep analyses to detect genetic variations associated with egg number based on a 55 K SNP Chicken chip in a population of 11,279 meat-type chickens.

Animals and phenotypes
A total of 11,279 meat-type chickens were obtained from seven lines of two breeding companies in China. The W1, W2, W3, and W4 were four commercial white-feathered pure lines in Yunnan Province, China. Body weight exceeded 2.4 kg at 40 days of age in lines W1, W2, and W4. Body weight exceeded 1.6 kg at 40 days of age in line W3. The feed conversion ratio was under 1.69 between 28 and 40 days of age, and the age at first egg was 175-180 d.
The Y1, Y2, and Y3 were three commercial yellow-feathered pure lines in Jiangsu Province, China. The feed conversion ratio was ~ 3.13, which was considered medium growth rate meat-type chickens.
Line W2 was selected for the egg production more than eight generations. We used the individuals of three generations (3159 chickens: 2532 females; 627 males), including the 6th, 7th and 8th generations. All other lines were used for individual data of one generation. All chickens had pedigree, genomic information and were housed in identical individual cages. There were 11,279 individuals in all lines, of which 9279 were recorded for egg numbers. Blood samples were collected from the brachial veins of each chicken using the standard procedure of the breeding program, which was approved by the Animal Welfare and Ethics Committee of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (IAS-CAAS, Beijing, China).

Staged egg numbers statistics
We made a detailed division of the laying period using 627 hens from the 7th generation of line W2. We recorded the age at first egg (AFE), and egg production for each hen from first egg to 354 days of age. We divided the egg-laying period into four stages according to the laying curve in Additional file 1: Fig. S1 and counted the egg numbers (EN) at each stage. Four stages, including the stage of rapidly increasing egg production (from the onset of laying eggs to 195 days of age; EN1); the stage of peak egg production (from 196 to 227 days of age; EN2); the stage at which the laying rate was between 60% and 75% (from 228 to 307 days of age; EN3); and the stage of decreasing egg production (from 308 to 354 days of age; EN4). Moreover, we counted the total egg number from first egg to 354 d (TEN).

Genotyping and quality control
Genomic DNA was extracted from blood samples using the phenol-chloroform method. All the birds were genotyped using a customized 55 K SNP chicken array obtained from Beijing Compass Biotechnology Co., Ltd. (Beijing, China) [9]. We updated the SNP positions using the newest release from the University of California-Santa Cruz (GRCg6a/galGal6 genome version). The following quality control criteria were applied to the target panel: individual call rate ≥ 90%, SNP call rate ≥ 90%, and minor allele frequency ≥ 0.05. Ultimately, 37,847-45,346 SNPs and 11,279 birds remained for further analyses.

Population structure testing
Twenty individuals randomly selected from each line for structure using ADMIXTURE software (version1.3.0) [10]. A principal component analysis (PCA) plot was constructed using the "scatterplot3d" package in R (version 4.0.2).

Genome-wide association study
GWAS was performed using 9778 individuals with egg number phenotypes from the all lines. The GWAS for egg number traits was performed using the univariate linear mixed model implemented in GEMMA software, version 0.98.1 [11]. The statistical model was as follows: where y is the vector of phenotypic values; W is the vector of covariates, including a column of 1 s; α is the vector of the corresponding coefficients, including the intercept; x is the vector of marker genotypes; β is the effect size of the marker; u is the vector of random polygenic effects; ε is the vector of errors; τ −1 is the variance of the residual errors; λ is the ratio between the two variance components; K is the centered relatedness matrix estimated from SNPs, and Ι n is the identity matrix. MVN n is the n-dimensional multivariate normal distribution. The Wald test was used to select SNPs associated with egg production traits. We used the parameters of plink software --indeppairwise 25 5 0.2 to infer effective independent tests. Considering the over-conservation of the 5% Bonferroni correction method, we adjusted the threshold of genome-wide significant P-values based on the number of independent SNPs [12]. Manhattan and Q-Q plots were constructed for each trait using the "qqman" package in R (version 4.0.2). Boxplots were produced by the "ggplot2" package in R (version 4.0.2).

Meta-analyses
Meta-analyses of the seven lines were performed using standard method by Metal software (version 2011-03-25). Meta-analyses can combine either (a) test statistics and standard errors or (b) P-values across multiple GWAS for a single trait (taking sample size and direction of effect into account) [13]. The detailed calculation formula is introduced in reference [13]. Meta-analysis results in little or no loss of efficiency compared to analysis of a combined dataset including data from all individual studies.

Association study on GGA Z
Due to the particularity of GGA Z, we used "-method threshold" function of SNPTEST v2.5.6 software [14], a software for sex chromosome association study, for GGA Z association study. For chickens, female genotypes are (1) y =W + x + u + ; u ∼ MVN n 0, −1 K , ∼ MVN n 0, −1 I n encoded as 0/1 in SNPTAST. In addition to association test statistics, SNPTEST can output expected genotype and allele counts for diploid and haploid samples. Computation of allele frequencies and info statistics also take into account ploidy.

Selective sweep analysis
Selective sweep analysis was performed on all hens from the line W1 and 782 hens from the 8th generation of line W2, which had genomic information obtained by 55 K SNP chips. Line W2 was selected for the egg number more than eight generations, Line W1 was not selected for egg number. We used the "--weir-fst-pop" function of the Vcftools v0.1.14 software to calculate the population differentiation index (Fst), which were a window length of 40 kb and a step size of 2 kb [15].

Bayesian analysis
Two Bayesian approaches, Bayes B and Bayes R, were used to obtain the estimated marker effect explained by each SNP. Bayes B and Bayes R analyses were performed using the BGLR and hibayes package in R, respectively. The number of iterations after the burn-in phase and that of the burn-in period were 20,000 and 10,000, respectively. The purpose of the Bayesian analysis was to estimate whether the QTL we mapped had significantly larger effects. The Bayesian analysis' result was the effect size of the SNPs, while the result of GWAS was the P value and effect size of SNPs associated with phenotype. Their results could be verified by each other.

Genome prediction models
The genome prediction models, genomic best line prediction (GBLUP) and single-step GBLUP (ssGBLUP), were utilized in this study and were as follows: Where y is the vector of phenotypic value, X and Z are the association matrices of fixed and additive genetic effects, b is the vector of the fixed effect, α is the vector of the random additive genetic effect, and e is the vector of the random residual error. For GBLUP, when an α obeyed the following normal distribution α ~N(0, G σ 2 u ), G was the consanguinity matrix (G matrix). In ssGBLUP, the H matrix replaces the G matrix.
H matrix construction: Where H is an inter-individual relational matrix using genome-wide markers and pedigrees, A is the pedigree-based relatedness matrix, G is the genomebased relationship matrix, and A 22 is the individual genealogical relation matrix of sequencing individuals. The relative weights of G and A 22 in the H matrix were set as The level of the G matrix was corrected to the A 22 (matrix) G * represents the adjusted G matrix, and a and b are as follows: The formula of the H matrix after merging is: The animal models of the two traits are as follows: Where y, X, b, Z, α, and e are the same formula (2). The genetic parameters of line W2 for egg production traits were calculated using the ssGBLUP model, and the heritability of the seven lines was estimated by the GBLUP model. GBLUP and ssGBLUP were calculated using ASReml v4.1 software [16].

Basic statistics and genetic parameters
A total of 9778 hens from seven lines were subjected to genome-wide association analyses. The egg number of each line was evaluated at a different time. The heritability of egg number varied between the seven lines. The results showed that the egg number traits exhibited low heritability (0.034-0.258) ( Table 1). The basic statistics for the distribution of egg production traits from the 627 hens used in the GWAS and genetic parameters are shown in Table 2. EN1 and EN4 had the largest phenotypic variation, as some chickens were late to laying and laid eggs at inconsistent intervals. Additionally, some birds died or were too old to lay eggs in EN4 because of their low rate of follicle development   [17]. A total of 7771 independent tests of all chromosomal SNPs were performed. Therefore, the threshold P-values for genome-wide suggestive and genome-wide significance were calculated as 1.29E-04 (1.00/7771) and 6.43E-06 (0.05/7771), respectively. The heritability and correlation of egg production traits was estimated using a relationship matrix composed of pedigrees and genotypes. AFE and EN1 had medium heritability (0.24 and 0.27), while the heritability of ENs in other egg-laying stages was low, ranging from 0.09 to 0.17. AFE had a high negative correlation with EN1, a moderately negative correlation with TEN, and a low correlation with ENs in other egg-laying stages. It was shown that AFE affected the laying performance by affecting the EN1. EN1 had a moderately positive correlation with TEN, and low positive correlations with ENs in other egg-laying stages. EN2, EN3, EN4, and TEN exhibited high correlations with each other (> 0.8) ( Table 3).

Population structure testing
PCA using the first three principal components showed obvious population stratification between the seven lines and that four lines (W1, W2, W3, and W4) of whitefeathered meat-type chickens were mixed with each other (Fig. 1). Figure S2 displays a bar plot based of the cross-validation error rate. When K = 2, the fast-growing whitefeathered and the yellow-feathered meat-type chickens appeared as two differentiated clusters. When K = 3, the four white-feathered meat-type chicken lines were divided into three groups. When K = 4, the three yellowfeathered populations were separated (Additional file 3: Fig. S3).

GWAS for single-line egg number traits
GWAS was performed on each line, and 35 SNPs were associated with egg number. Using the Ensembl database to annotate SNPs, we found 31 genes around significant SNPs (Additional file 6: Table S1). For the GWAS of the seven lines, the genomic inflation factors (λ, 0.993-1.069) were calculated to evaluate the accuracy and reliability of the GWAS results (Additional file 4: Fig. S4).

Meta-analysis based on GWAS results of single-line egg number traits
Through the above analysis, the SNP effect of egg numbers in a single line could be obtained, but it was unknown whether the same region exists for multiple lines with large genetic background differences. Therefore, using the above GWAS results, a meta-analysis was performed on all individuals of the seven lines, and an obvious bulge region was found on the Z chromosome ( Fig. 2A).

Significant region was differentiated between lines with different breeding directions
The Fst and major allele frequency analyses were performed on the line W1 and W2 chickens ( Fig. 2B and C). It was a significant signal of selection that a QTL region overlapped with the result of meta-analysis located between 10.81 Mb and 13.05 Mb on the GGA Z. The region contains 36 genes. Eight genes were identified as candidate genes after functional annotation. The genes included glial cell derived neurotrophic factor (GDNF), DAB adaptor protein 2 (DAB2), protein kinase AMPactivated catalytic subunit alpha 1 (PRKAA1), NAD kinase 2, mitochondrial (NADK2), WD repeat domain 70 (WDR70), leukemia inhibitory factor receptor alpha (LIFR), complement C6 (C6), and complement C7 (C7).

Phased egg number traits to verify significant QTL region
(1) GWAS for staged egg number traits Combined with the results of GWAS and selective sweep analysis, we further analyzed the region between 10.81 Mb and 13.05 Mb on the GGA Z in line W2. For the GWAS of 627 hens in line W2, Manhattan and quantile-quantile (Q-Q) plots for egg production traits were generated, as shown in Fig. 3 and Additional file 5: Fig. S5. Detailed information on the SNPs significantly associated with EN2, EN3, EN4, and TEN is shown in Table 4. For the EN2, 26 genome-wide significant SNPs were identified on GGA Z, and 38 genome-wide suggestive SNPs were detected on GGA Z, GGA1, GGA3, GGA4, GGA10, and GGA13. The genomic inflation factor (λ) was 1.00 for EN2, which suggested that the population structure was well controlled.
For the TEN, one significant SNP and 38 suggestive SNPs were detected on GGA Z, GGA1, GGA3, and GGA23. The SNPs on the GGA Z were also associated with EN2. Three SNPs (rs318154184, rs13769886, rs313325646) associated with EN2, EN3 and TEN, were located on or near the prolactin receptor (PRLR) gene on the GGA Z. There were 33 coincidently associated SNPs for EN2 and TEN between 10.81 Mb and 13.05 Mb on GGA Z, among which 29 SNPs were associated with significantly different egg number EBVs of the 3159 individuals in line W2. The 29 SNPs were located in four LD blocks (Fig. 4), and the visualization for analysis of variance of one SNP selected each Block.
(2) Association tests on GGA Z Considering the particularity of GGA Z, we used SNPTEST software to analyze GGA Z of EN2 and TEN traits of 627 hens in line W2 specially. The results (Fig. 5) were consistent with those in (1), which could indicate the reliability of the region between 10.81 Mb and 13.05 Mb on GGA Z.

(3) The effect of each SNP
We used Bayes B and Bayes R to estimate the effect value of each SNP using the EN2 and TEN phenotypes of 627 hens from line W2. There were obvious bulges in the significant QTL region, which was similar to GEMMA in that the calculated effect size also had obvious bulges in the region (Additional file 6: Fig. S6). In the effect estimated for this region, GEMMA, Bayes B and Bayes R estimated EN2 and TEN accounted for 0.89%, 1.01%, 0.47%, 0.37%, 0.52%, and 0.52% of total effects, respectively. This result indicated that there were indeed strong effects in the region that could affect phenotypic variation.

Significance analysis between phenotypes and SNPs of the QTL region
To further verify the correlation between egg numbers and 29 SNPs, we analyzed the SNPs' change of dominant genotypes with breeding. We found that the frequency of dominant genotypes gradually increased with the selection of W2 for egg production, which indicated that the SNPs were also selected with the selection of phenotype, which could well explain these SNPs and egg production is significantly correlated (Fig. 6). Therefore, the 29 SNPs were identified as important candidate SNPs (Table 5).

Discussion
In the twentieth century, layers and meat-type chickens were developed in different directions to avoid conflicts between egg-producing and meat-producing traits [18]. Such activities were extremely successful in improving productivity; laying hens produced more than 320 eggs during their 52-week laying period, while broilers gained 50-60 times their weight from hatching to sale [19]. Those activities also led to substantial genomic differences between layers and meat-type chickens.
To explore the genetic mechanisms of laying traits in meat-type chickens, 11,279 chickens were used in this experiment.
The genetic contributions of egg number have low or medium levels of heritability, typically ranging from 0.16 to 0.64 [6]. However, among the seven lines evaluated in this study, the heritability of egg number was different (0.034-0.258), but it was not the color of their feathers that caused the difference. The range of heritability observed in this study was consistent with those reported in White Leghorn hens (0.05-0.27) [20]. The main reasons for other differences in heritability included environmental effects, population size, animal breeds, differences in assessment methods, and the reliability of the individual information records.
Further analysis the EN of partial periods and the total period in Line W2, Liu et al. predicted heritability range between 0.05-0.24 for White Leghorn hens [12]. Our results showed that AFE and EN1 had a high level (0.24-0.27) of heritability, with the heritability of the other stages being 0.093-0.17, which was slightly lower than those previously reported (0.17-0.20) [21]. AFE had a strong negative genetic correlation with EN1 (− 0.96), and a moderate negative correlation with TEN (− 0.37), which was consistent with the findings of Yuan et al. [21] and Liu et al. [12].
Meta-analyses have been widely used in multi-population GWAS studies. Falker-Gieske et al. used 1306 White Leghorn chickens and a meta-analysis of GWASs to identify the variations and genes related to feather pecking, and found 15 significant SNPs that may be related to feather pecking [22].
We used GWAS, meta-analyses, selective sweep analyses, and Bayesian analyses to scan the SNPs and QTLs associated with laying traits in meat-type chickens. A 10.81-13.05 Mb region on GGA Z was analyzed using these methods and 36 genes were annotated in this region. However, no studies have detailed the genes related to egg number traits in chickens. Notably, several regions associated with egg number traits have been identified on the GGA Z (https:// www. anima lgeno me. org/ cgi-bin/ QTLdb/ GG/ index). The regions identified in our study were 1.9 Mb apart from those reported by Zhao et al., which were carried out in more than a dozen breeds, including White Leghorn and Rhode Island Red chickens [23]. Although the results were different between the present study and that of Zhao et al., they suggest that GGA Z is an important genetic region involved in egg production.
Referring to prior studies, we found that the regions between 10.81 Mb and 13.05 Mb on the GGA Z identified in this study were inconsistent with those reported in layers, which may be due to the differences in the genomes between layers and meat-type chickens. In the regions, genes that may be involved in reproductive traits were identified, including the GDNF gene, which has been reported to affect the in vitro maturation of human oocytes, increasing the prevalence of MII stage oocytes from 34% to 49% [24]. GDNF has also been reported to affect the development of oocytes in a variety of animals, such as humans, pigs and mice [25][26][27].
DAB2 is an endocytic adaptor protein in several NPXY motif-containing cell surface receptors, including the lipoprotein receptor and integrin. Budna et al. also suggested that the DAB2 gene may be involve in the regulation of MII phase oocyte formation and other processes critical to porcine fertility [28]. The PRKAA1 gene encodes the catalytic α-subunit of 5′ AMP-activated protein kinase has an opposite effect--maintenance of the meiotic block in porcine and bovine oocytes [29,30] AMPK can control steroidogenesis in ovarian cells (granulosa, theca cells and corpus luteum cells) and germ cell maturation [31]. The NADK2 gene in this region, the NAD kinase (NADK) is the sole NADP biosynthetic enzyme [32]. NADPH, which is the reduced form of NADP, regenerates cellular oxidative defense systems via glutathione reductase to counteract oxidative damage [33,34]. Therefore, NADK2 may affect egg production by affecting reactive oxygen species levels [34]. The WDR70 gene was identified as a candidate gene for milk production and reproductive traits in Chinese and Northern European Holstein dairy cows [35]. The gene is homologous in chickens and cattle and may have the same function [36].
The LIFR gene can form polymeric complexes with other receptors, such as glycoprotein 130 (GP130), to stimulate the JAK/STAT, MAPK and PI3K signaling pathways [37]. JAK/STAT and PI3K signaling pathways influence the hypothalamic-pituitary-ovarian (HPO) axis to regulate the breeding cycle of laying hens. Two genes, C6 and C7, were identified in this region. The complement system plays an important role in the process of inflammation and infection, and also provides a link between innate and adaptive immunities. Immune capacity can affect the reproductive performance, health, and welfare of hens [38]. In addition, three SNPs (rs318154184, rs13769886, rs313325646) associated with EN2, EN3, and TEN were located on or near the prolactin receptor (PRLR) gene on GGA Z. The PRLR gene is associated with laying traits and consists of 15 exons and 14 introns in chicken [39,40]. Liu et al. hypothesized that PRLR might be the primary gene responsible for sexual maturation and is considered a candidate genetic marker for reproductive traits [41]. Studies showed that the prolactin gene is expressed in the hypothalamus, pituitary, oviduct, and ovary of chickens. Some studies also showed that polymorphisms occur at different locations of the PRLR gene, including SNP in exons 2, 5 [40], and 10 [39], which were mainly associated with egg number traits.

Conclusion
In conclusion, our results indicated that the heritability of egg numbers was different in different chicken populations, and egg number was a low-heritability trait. Based on GWAS analyses of 627 individuals in one generation of line W2, meta-analyses of seven lines, and selective sweep analyses of lines W1 and W2, we found an interesting region on GGA Z. Where multiple SNPs and genes may be significantly related to egg production. Our findings contributed to a better understanding of the genetic basis of egg production.